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Abstract. We measured the radial-velocity curve of HD 77581, the B-supergiant companion of the X-ray pulsar Vela X-l, 
using 183 high-resolution optical spectra obtained in a nine-month campaign. We derive radial-velocity amplitudes for different 
lines and wavelength regions, and find all are consistent with each other, as well as with values found in previous analyses. We 
show that one apparent exception, an anomalously low value derived from ultra-violet spectra obtained with the International 
Ultraviolet Explorer, was due to an error in the analysis procedures. We re-analyse all IUE spectra, and combine the resulting 
velocities with the ones derived from the new optical spectra presented here, as well as those derived from optical spectra 
published earlier. As in previous analyses, the radial velocities show strong deviations from those expected for a pure Keplerian 
orbit, with root-mean-square amplitudes of ~ 7 km s _1 for strong lines of Si IV and N III near 4100 A, and up to ~ 20 km s _1 
for weaker lines of N II and Al III near 5700 A. The deviations likely are related to the pronounced line-profile variations seen 
in our spectra. Our hope was that the deviations would average out when a sufficient number of spectra were added together. 
It turns out, however, that systematic deviations as a function of orbital phase are present as well, at the 3 kms -1 level, with 
the largest deviations occurring near inferior conjunction of the neutron star and near the phase of maximum approaching 
velocity. While the former might be due to a photo-ionisation wake, for which we observe direct evidence in the profiles 
of US and Ha, the latter has no straightforward explanation. As a result, our best estimate of the radial-velocity amplitude, 
K op t = 21.7 ± 1.6 kms -1 , has an uncertainty not much reduced to that found in previous analyses, in which the influence of 
the systematic, phase-locked deviations had not been taken into account. Combining our velocity amplitude with the accurate 
orbital elements of the X-ray pulsar, we infer M nB sin 3 i = 1.78 ±0.15 M@***. 
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1. Introduction 

Neutron stars follow a mass-radius relation that depends on 
the equation of state of the superdense matter in their inte- 
riors. For densities up to the nuclear saturation density (~ 
10 14 - 5 gem -3 ), the equation of state is reasonably well con- 
strained by laboratory measurements, but for higher densities 
it is still under much debate. For a given equation of state, 
there is a maximum possible mass of a neutron star. For a 
"soft" one (i.e., "compressible" matter), the maximum is near 
1.5 Mq, while for a "stiff" one (i.e., incompressible matter) 
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it can be as high as 2 M Q (for recent reviews, see Heiselberg 



& Pandharipande 2000; Lattimer & Prakash 2000; Srinivasan 



2001). In this respect, the X-ray pulsar Vela X-l is of par- 



ticular interest, as it is the only neutron star for which there 
is strong evidence that it is more massive th an the c anonical 
1.4 Mq (for a re view, see Van Kerkwijk et al. 1995a ; Thorsett 
&Chakrabarthy|l999|). 



Vela X-l (4U 0900-40) is a 283-s pulsar orbiting the 
B0.5 lb supergiant HD 77581 in a period of 8.96 days. The 
X-ray source is powered by accretion from the B supergiant's 



wind (see Kaper 2001 for a recent overview of high-mass X-ray 



binaries hosting an OB-supergiant). Since the discovery of the 



X-ray pulsations in Vela X-l by McClintock et al. (1976), its 
slightly eccentric (e = 0.09) orbit has been determined with in- 
creasing accuracy from Doppler-shift measurements of pulse- 
arrival times. Recently, an update of the orbital parameters has 
been obtained based on the five-year long series of frequency 
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measurements in hard X-ray flux with the Burst And Transient 
Source Experiment (BATSE) on board the Compton Gamma- 
Ray Observatory (Bildsten et al. 1997); see Table |j. 

The evidence for VelaX-l's high mass is the relatively 
large amplitude, K opt , of the radial-velocity variations of the 
B supergiant companion, first noticed by Van Paradijs et al. 
(1976). From their analysis of close to a hundred photo- 



graphic spectra, Van Paradijs et al. (1977a, hereafter Paper I) 
found K opt = 21.75 ± 1.15 km"s~"\ 



which, combined 



with the BATSE orbit, corresponds to a neutron-star mass of 
M ns sin 3 i = 1.79 ± 0.11 M© (here M ns is the mass of the 
neutron star, and i is the orbital inclination; uncertainties are 
la unless specified otherwise). Guided by the prospect that 
a more accurate mass determination of Vela X-l might pro- 
vide an important constraint on the equation of state appro- 
priate for matter in neutron stars, Van Kerkwijk et al. ( 19951 , 
Paper II) made a detailed study of the radial-velocity curve of 
HD 77581, based on 40 high signal-to-noise CCD spectra, as 
well as 26 ultraviolet spectra obtained with the International 
Ultraviolet Explorer (IUE) and 13 nightly averages of digitized 
photographic spectra. They derived A" pt = 20.8±1.4kms -1 , 
confirming the value found in Paper I. 

In Paper II, the accuracy of K opt was found to be limited 
not by the statistical uncertainty in the radial-velocity measure- 
ments, but by systematic non-orbital contributions to the ra- 
dial velocity. These contributions are correlated within a given 
night and can reach amplitudes up to ~ 15 kms -1 . They did 
not appear to be correlated with orbital phase. It was suggested 
that these velocity deviations might be caused by non-radial 
pulsations of HD 77581, similar to those seen in many early- 
type stars, perhaps excited by the varying tidal force excerted 
by the neutron star in its eccentric orbit. The situation may be 
complicated further by systematic effects on the spectra due 
to a photo-io nisati on wake trailing the neutron star in its orbit 
(Kaperetal. [l994l ). 

In Paper II, it was suggested that for a better understand- 
ing of these radial-velocity excursions, and thereby an im- 
provement of the accuracy of the radial-velocity curve, either 
an extensive long-term sequence of spectroscopic observations 
would be required, with which one might "average out" the de- 
viations, or a continuous spectroscopy campaign with a time 
resolution of a few hours, extending over at least one orbital 
cycle. We followed the first approach and monitored the binary 
system for a period of about 9 months, taking high-resolution 
spectra of HD 77581 on a night-by-night basis at the European 
Southern Observatory. 

In Sect. ||, we describe the observations and the data- 
reduction procedures. In Sect. the spectra are presented. In 
Sects. |] and ||, we describe the method of determining ra- 
dial velocities and discuss the resulting radial-velocity curves. 
Sect. U is an intermezzo, dedicated to the reanalysis of IUE 
spectra, a subset of which was already used in Paper II in the 
determination of the mass; we compare our results to the ones 



found by Stickland et al. ( 1997 > who analysed the same set of 



IUE data, and explain the origin of the difference between our 
results. In Sect. [7j we combine all velocity curves. We also in- 
vestigate systematic deviations with orbital phase, finding that, 
unfortunately, these are present and that they limit our final ac- 



curacy. In Sect. ^| we present our results for the mass of the 
neutron star and in Sect. ^ we discuss how these results influ- 
ence the current status of our knowledge on the internal struc- 
ture of neutron stars. 

2. Data reduction and analysis 

We observed HD 77581, the bright (V = 6.9 mag) B- 
supergiant companion to Vela X-l, in the period October 15, 
1995, until July 15, 1996, using the 1.4-m Coude Auxiliary 
Telescope (CAT) and Coude Echelle Spectrograph (CES) at 
the European Southern Observatory (ESO), La Silla, Chile. 
The observations were carried out by the night assistant of 
the CAT (in La Silla or, more often, remotely from the ESO 
Headquarters in Garching, Germany) or by the astronomer 
present at the telescope. 

Every night one spectrum (in a few cases two) of HD 77581 
was obtained, as well as a few bias frames, two internal lamp 
flat fields and a Th-Ar wavelength calibration frame. Exposure 
times range from 20 to 30 minutes and the typical signal-to- 
noise ratio of the extracted spectra is il 300 per 5 kms -1 
resolution element. A fraction of the obtained spectra (47 out 
of 230) could not be used due to various practical problems 
(e.g., wavelength region covered by the stellar spectrum was 
not identical to that in the wavelength-calibration frame, no cal- 
ibration frame available, wrong exposure level science frame, 
etc.). 

The CES spectrograph has two "paths", one optimized for 
the blue and one for the red part of the optical spectrum. The 
CES long camera is equiped with one of two CCD detectors, 
CCD#34 (2048 x 2048 pixels) orCCD#38 (2720 x 520 pixels, 
and thus a larger wavelength coverage). In a given night, one 
has to select the red or the blue path, as well as the CCD. In or- 
der to minimize the impact of our programme on the observing 
run carried out during the remaining part of the night, we had to 
observe with the setting selected by the observer of that night. 
For this reason, we selected two wavelength regions, one in the 
blue (4077^1-127 A) and one in the red path of the spectrograph 
(5670-5720 A). A total of 104 usable spectra were taken in the 
red and 79 in the blue wave leng th re gion. A log of the observa- 
tions is given in Tables A 1 and A2 in the Appendix. 

One of the main problems we were faced with during the 
analysis is the vignetting in the science frames, which varies 
from night to night, depending on the position of the telescope. 
It is caused by the small size of the field lens of the long cam- 
era (which was not designed for large CCDs), and the varia- 
tions are due to the large distance between spectrograph and 
telescope. The CAT telescope is positioned in a dome just next 
to the dome of the 3.6m telescope, while the CES spectrograph 
is located inside the 3.6m dome. The light coming from one of 
the Nasmyth foci of the CAT goes through a pipe connecting 
the two domes, into the predisperser area of the spectrograph, 
bridging a distance of about 15 m. For this reason, it is very dif- 
ficult (especially in view of the regular earthquakes) to align the 
telescope and spectrograph beam. Dome flats provide a better 
way to correct for this misalignment than internal flats; unfor- 
tunately, in practice it was not feasible to obtain dome flats for 
our programme. The vignetting is more severe at longer wave- 
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Fig. 1. The average blue CAT/CES spectrum centered at HS 
4102 A, with line identifications. The wavelength scale is in 
An gstr0m. 



Fig. 2. The average red CAT/CES spectrum centered at 5695 A, 
which includes several lines of a N II multiplet. Note the differ- 
ence in strength between the lines present in the blue and red 
spectra. 



lengths. The correction for the vignetting, though hampered by 
the relatively small wavelength region covered in one spectrum, 
has been carried out during the normalisation procedure. 

Below, we describe the extraction, calibration, and normal- 
isation of the spectra. 

2. 1. Bias subtraction 

The electronic bias distribution is uniform over both chips, 
therefore in the bias correction (per night) we subtracted a num- 
ber, derived either from the mean of the bias frames taken in 
that night (for the spectra taken from October 95 until March 
96), or from the overscan of the science frame itself (for the 
spectra taken from March 96 until July 96). 

2.2. Flatfielding 

In order to apply the flatfield correction, we constructed one 
"master" flatfield frame for each of the two CCDs and the 
two wavelength regions, i.e., four frames in total. The mas- 
ter flatfield consists of the average of all available flatfields in 
a given wavelength region and taken with a given CCD, after 
inspection for saturation or under-exposure and a consistency 
check of their spectral parameters (central wavelength, resolu- 
tion power, pre-disperser exit slit width). The individual flat- 
fields were normalised to unity before averaging. Since the po- 
sition of the spectrum on the CCD was not always the same dur- 
ing the extended period of our observing campaign, the master 
flatfield does not have a uniform exposure level. We scaled the 
flatfield to unity by dividing each row of the detector by its 
mean value. The night-to-night changes in instrument setting 
are small and have no significant effect on the flatfield struc- 
ture (e.g., we were not confronted with strong detector fringes 
or ghost images). 



2.3. Optimal extraction 

We used an op timal extraction algorithm devel oped by Van 
Kerkwijk ( 1993 ), based on the method of Home ( 1986 ), which 
includes corrections for small non-alignments of the spectrum 
with the dispersion direction of the chip, and detects and fil- 
ters for cosmic-ray events. Some bad columns on the detectors 
have been corrected for by interpolation (in pixel space) over 
the gaps in the extracted spectrum. 

In the red wavelength region some atmospheric water 
vapour lines are present. As they would cause problems in the 
velocity determination we removed them by linear interpola- 
tion. 



2.4. Wavelength calibration 

For the wavelength calibration we use the Th-Ar calibration 
frames, fitting the line positions with a second-order dispersion 
relation. The typical root-mean-square residuals are 0.001 A. 
The dispersion is 0.018 A per pixel for the blue region and 
0.024 A per pixel for the red region. There are about five pixels 
per resolution element. 



2.5. Barycentric correction and orbital ephemeris 

We corrected for the Earth's orbital and rotational motion by 
calculating the component of the Earth's velocity in the di- 
rection of HD 77581 at the mid-time of the observation and 
subtracting it from the velocity measured from the determined 
Doppler shift of the line (see Sect. 4). The orbital phase was 
calculated using th e orbi tal parameters based on the BATSE 
data (Bildsten et al. |1997[ ); these are listed in Table [| 



2.6. Normalisation of the spectra 

Normalising the spectra was not a straightforward procedure, 
the main reason for this being the (time-dependent) vignetting 
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Table 1. The wavelength regions used to normalise the spectra, 
for the blue and red wavelength region and CCD#34 and 38, 
respectively. Before normalisation, the selected continuum re- 
gions are shifted by an amount corresponding to the (predicted) 
Doppler shift for the orbital phase at the time of observation. 



CCD#38 



CCD#34 



Blue 












4080 


45-4082 


10 


4086 


04-4086 


50 


4091 


25-4091 


77 


4091 


50-4092 


50 


4107 


10-4108 


00 


4094 


35-4094 


95 


4122 


62^1124 


05 


4108 


10-4108 


50 








4113 


20-4113 


70 








4118 


22-4118 


77 


Red 












5672 


30-5672 


90 


5673 


40-5673 


55 


5682 


82-5683 


38 


5682 


82-5683 


38 


5693 


00-5693 


66 


5693 


00-5693 


66 


5699 


98-5702 


58 


5699 


61-5700 


38 


5714 


25-5716 


00 


5713 


04-5714 


04 



(see above) which affects all spectra. We explored different 
methods to normalise the blue and the red spectra, taking into 
account: (i) the use of the same continuum regions (after shift- 
ing the spectra approximately to the reference frame of the star) 
for each spectrum and (ii) the assumption that the shape of the 
spectral continuum can be described by a polynomial of order 
four or lower. 

The main problem we were confronted with was the lack 
of clean continuum regions due to the small wavelength region 
(~ 50 A) covered in each spectrum. Also, the intrinsic (i.e. non- 
orbital) variability of the stellar spectrum sometimes causes the 
lines to extend into the spectral continuum regions. 

Below we describe the normalisation methods we preferred 
after testing different approaches. We note that we repeated the 
full analysis (including the cross-correlation procedure) for dif- 
ferent normalisation methods, but did not find significantly dif- 
ferent results. 

2.6.1. Blue spectra 

The central part of the blue spectrum (Fig. [j]), i.e., at the blend 
of N III 4097 A and H<5, is not affected by vignetting, so the 
shape of the continuum is well reproduced by a linear function 
defined by the continuum regions left and right of the blend. 
The edges of the spectrum were normalised using two second- 
order polynomials, one at each side of the spectrum, which 
were forced to have the same slope as the linear function at 
the point where they meet. The normalisation regions are given 
in Table | 

2.6.2. Red spectra 

The effect of vignetting is more pronounced in the red wave- 
length region. Furthermore, the intrinsic variability of the line 
profiles is much larger compared to that in the blue spectrum, 
the spectral lines are closer to each other, and sometimes there 



appear to be (variable) P Cygni emission features at the blue 
side of the strongest N II lines. There are no extended contin- 
uum regions which are completely devoid from spectral fea- 
tures. The regions listed in Table [j] are those we believe are 
least affected by spectral lines. 

In order to reduce the effect of the vignetting (which other- 
wise would force us to use a high-order polynomial to fit rea- 
sonably the shape of the continuum) we divided all the spec- 
tra by one spectrum (r8005r) which appears to be only mod- 
erately affected by vignetting. The residuals were fitted with 
a third-order polynomial. The latter fits should approximately 
reproduce the shape of the continuum and eliminate most of 
the vignetting. Then, we multiplied the original spectra with 
this fit. In this way, the spectra were "ironed". Finally, we fit- 
ted a fourth-order polynomial to the ironed spectra, using the 
normalisation regions listed in Table |} 

3. Description of the spectra 

In Figs. |l] and || we show the average of the blue and of the red 
spectra, respectively, after having taken into account the dif- 
ference in orbital motion. The blue spectrum includes stronger 
lines than the red spectrum, which makes the blue spectrum 
better suited to measure the star's radial velocity. 

An overview of the blue and red spectra is presented in 
Figs. |]and |[ respectively. The spectra are ordered according to 
orbital ph ase; the c orre sponding observing times can be found 
in Tables |A.l| and A2 in the Appendix. The few spectra ob- 
tained with CCD#34 (instead of CCD#38) are shown in Figs. | 
and (7|. The wavelength scale of the spectra is with respect to the 
stellar rest frame, i.e., the orbital motion has been taken out to 
facilitate comparison. The figures show that the profiles of the 
lines in the red spectra are much more variable than those of the 
lines in the blue spectra. In a given spectrum, the different lines 
seem to vary mostly in concert. In some cases the variations are 
dramatic, especially in the red spectra, where the depths of the 
photospheric lines are sometimes strongly reduced. The vari- 
ations do not show an obvious dependence on orbital phase. 
In order to have a better view on possible orbital-phase related 
line-profile variations, we binned the spectra in nine phase bins 
(nine bins is chosen since the orbital period is very close to nine 
days). The results are shown in Figs. || and ^| 

3. 1. The blue spectra centered on H5 4102 A 

The blue spectra cover HS and two strong Si IV lines, at 
4088.85 and 41 16. 10 A, as well as weaker lines of N III (at 
4097.31 A), O II (41 12 A) and He I (4121 A). Furthermore, 
there are two small features which we could not identify; these 
lines are also seen in the average spectrum in Paper II. These 
are located at 4084.88 A just shortward of the Si IV 4088.85 A 
line and at 4092.96 A at the opposite side of the Si IV line. 

In general, the lines show strikingly similar behaviour. For 
the two Si IV lines this is expected, since they arise from the 
same multiplet. But the N III fine at 4097.31 A, while blended 
with the HS line, generally behaves similarly as well. See, e.g., 
the blue spectrum at phase 0.322 in Fig. |l where the profiles of 
the Si IV lines and the N III line show an identical deformation 
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Fig. 4. The blue spectra obtained with CCD#34, ordered ac- 
cording to orbital phase. 




Fig. 5. Mean blue spectra, obtained after averaging over nine 
different phase bins. The dotted line represents the mean profile 
at orbital phase 0.01-0. 1 1 and is plotted for comparison. 

of line shape. It is less easy to judge the case for the O II line, 
which is rather weak, and the He I line, which is not completely 
covered by our spectra (hence we did not include the latter in 
our analysis). 

Also in the phase-binned spectra, the lines appear mostly 
to vary in concert. The profiles close to X-ray eclipse (in phase 
intervals 0.01-0.1 1 and 0.95-1.00) are rather symmetric, while 
the mean profiles in phase intervals 0.17-0.33 and 0.75-0.80 
are shifted to the red, and those from 0.53-0.64 appear asym- 
metric to the blue. The variations are quite a bit smaller, how- 
ever, than the large internal spread within each phase interval. 

The only line showing deviant behaviour is HS. During or- 
bital phases 0.45-0.65, i.e., when the neutron star passes in 
front of the supergiant companion, the line profile is strongly 
distorted, especially the blue wing, where often a weak and 
blue-shifted absorption component can be observed, which can 
have such a high, blue-shifted velocity that it blends with the 
N III line. At other phases, the line profile behaviour is more 
similar to what is seen for the Si IV and N III lines. 

In Fig. ||, the enhanced blue absorption in H<5 can be seen 
clearly in the mean profiles of phase intervals 0.40-0.46 and 
0.53-0.57. This distortion of the HS line is very similar to that 




Fig. 7. The red spectra obtained with CCD#34, ordered accord- 
ing to orbital phase. 




Fig. 8. Mean red spectra, obtained after averaging over nine dif- 
ferent phase bins. The dotted line represents the mean profile 
at orbital phase 0.01-0.1 1 and is plotted for comparison. 

see n in t he H/3 and He I 4471 A line of HD 77581 (Kaper et 
al. 1994), and is most likely caused by the perturbation of the 



B-supergiant wind due to the presence of the X-ray source (see 



Sect. 7.2). 



3.2. The red spectra centered at 5695 A 

The red spectra cover lines which are much weaker than those 
in the blue, and which, as is clear from the spectra, also show 
much stronger intrinsic perturbations of their profiles. The lines 
covered include N II lines at 5676.02, 5679.56, and 5686.21 A, 
an Al III line at 5696.47 A, and an C III line at 57 10.77 A. 
Furthermore, the spectra contain a diffuse interstellar band at 
5705 A. Its constant profile demonstrates the stability of the 
instrument, but at the same time, since it is contaminating the 
C III line it hinders us from using the latter line to measure 
velocities. Finally, between the N II and the Al III line, at 
A569 1 . 1 A, there is a very shallow variable absorption feature 
which we are not able to identify. 

The N II lines belong to the same multiplet. All three lines 
show dramatic intrinsic variations, with, like for the blue spec- 
tra, little evidence for a correlation with orbital phase. On oc- 
casion, the distinction between the two blended N II lines dis- 
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Fig. 3. The blue spectra obtained with CCD#38 ordered according to orbital phase. The Doppler shift due to the orbital motion 
has been taken out. Note the variability of the line profiles, and the strong correlation in shape between the different lines. The H5 
line shows stronger variations in the blue wing when the X-ray source passes through the line of sight (0 = 0.5). These variations 
are likely due to the presence of a photo-ionization wake in the B-supergiant's wind. 



Fig. 6. The red spectra (CCD#38) ordered according to orbital phase. The lines of the N II multiplet show strong intrinsic 
variability, apparently unrelated to orbital phase. In the most extreme cases, the lines almost disappear. 
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appears completely (see, e.g., phases 0.184 and 0.511), or an 
additional feature occurs (phases 0.337 or 0.458), in an appar- 
ently irregular manner. The strongest N II line sometimes ap- 
pears to develop a red (P-Cygni) emission bump. If real, this 
would be an indication that the strongest N II line is formed 
in the outer expanding layers of the atmosphere. Note that the 
N II lines are often a useful diagnostic in studying non-radial 
pulsations in B-type stars (D.R. Gies, private communication). 
Fig. ^indicates that, on average, at orbital phases close to X-ray 
eclipse (0.95-1, 0-0.25) the red spectral lines are deepest. 

In general, the Al III and C III line profiles are similar to 
those of the N II lines. See, e.g., the spectra at orbital phases 
0.868 and 0.872, where all the line profiles (except that of 
the interstellar band, of course) are distorted towards longer 
and shorter wavelengths, respectively. Again, at several orbital 
phases (e.g., 0.01) features appear which could be interpreted 
as being due to incipient emission. These can be present ei- 
ther on the left wing (Al III line, e.g., phases 0.01, 0.296, and 
0.675), or on the right wing (Al III line, phase 0.448). In view 
of the difficulty of defining the continuum level, it is hard to say 
whether these are indeed emission features or whether there is 
a weak absorption feature right next to it. 

Sometimes (see, e.g., the spectra obtained at phases 0.653 
and 0.78) the depth of the lines decreases so much that the 
N II blend becomes very shallow, whereas the weaker N II 
and Al III lines almost disappear. 

4. Radial-velocity determination 

We determined the radial velocity of HD 77581 by means of 
cross-correlation. We used the same algorithm as in Paper II 



Table 2. Cross-correlation regions. 



CCD#38 



CCD#34 



(see also Van Kerkwijk 1993) where each spectrum is cross- 



correlated with all other spectra in the data set. We decided to 
use all the available spectra as independent reference spectra, 
because of the lack of a template spectrum of a bright star with 
identical spectral type. We did not cross-correlate against the 
average spectrum to avoid auto-correlation effects which lead 
to a systematic decrease of the amplitude of the radial-velocity 
curve. For details on the cross-correlation method we refer to 



Paper II and Van Kerkwijk (|199J). 

We first transformed the spectra to a logarithmic wave- 
length scale so that a Doppler shift is a linear displacement 
along the spectrum. The top part of the cross-correlation func- 
tion (CCF) was fitted with an analytic function (a Gaussian 
plus a linear function) whose maximum determines the velocity 
shift of the spectrum. Obviously, the relevant part of the cross- 
correlation function is the area close to the maximum; we fit 
only the upper 30 to 40% of the CCF, depending on its shape 
and smoothness (Fig. ^|). 

We cross-correlated all the lines separately as well as com- 
bined, using the wavelength regions listed in Table 0. The mea- 



sured radial velocities are listed in Table A. 1 (blue wavelength 



region) and Table A. 2 (red wavelength region) in the Appendix. 
The typical intrinsic, formal accuracies are 0.7 kms" 1 for the 
blue spectra and 1.5 kms -1 for the red spectra. 

The cross-correlation method worked satisfactorily for the 
blue spectra. The CCF is very smooth and the fit well de- 
fined. This results in an accurate determination of the veloci- 



Blue 

4084.40-4092.45 
4092.50-4099.04 
4108.30-4118.12 
Red 

5669.76-5682.98 
5682.98-5693.78 
5693.78-5701.91 



4085.00-4093.30 
4093.30-4099.40 
4110.13-4119.12 

5673.20-5682.98 
5682.98-5693.78 
5693.78-5701.91 



ties. In Fig. ^, we show a typical CCF, derived from the cross- 
correlation of the Si IV line, as well as the fit. 

In the red spectra, the CCF is much less well-defined, since 
the absorption features are much weaker. We partly overcame 
this by allowing the fit to the CCF to cover a larger area around 
the peak so that the determination of the peak position is not 
too much affected by the noise. Fortunately, the N II blend (at 
5676.02 and 5679.56 A), despite its strong intrinsic variations, 
is sufficiently deep to produce a reasonable CCF. In Fig. [j[ we 
show a typical CCF of the N II blend and of the Al III line. 

The presence of the interstellar line at 5705 A allows us 
to check for trends, e.g., due to instrumental effects, in the 
derived radial velocities. In the bottom-right panel in Fig. ^ 
we show the CCF for the diffuse interstellar band. It is rather 
wiggly, since the band is quite weak. Fig. |l0| displays the ra- 
dial velocities derived from the band (filled circles, displaced 
by —75 km s _1 ). For comparison, the radial-velocity curves of 
the N II 5686 A and Al III line regions are shown. The ra- 
dial velocity of the DIB is consistent with a constant value 
(Xicd — 1-1). A slight phase dependence might be introduced 
due to the neighbouring C III line, which moves close to the 
DIB around phase 0.75, but this does not become apparent in 
Fig. [l(^. Nevertheless, to be on the safe side, we did not use the 
radial velocities measured from the C III line in our solutions. 

From the above, we conclude that our error estimates are 
reasonable for the interstellar band. This feature does not vary 
in shape, and hence it satisfies the assumption underlying cross- 
correlation, viz., that the two spectra being correlated differ 
only in their radial velocity. This assumption will not hold 
for the stellar lines, given the large variations seen in the line 
profiles. Therefore, the error estimates become somewhat ill- 
defined. Nevertheless, if all lines varied truly in concert, one 
would expect the inferred velocities to be the same within the 
errors. 

In Fig. [llL we compare the radial velocities derived for the 
different lines. For the two regions centred on the Si IV lines, 
we find a good correspondence. However, while there are no 
systematic deviations, the velocities inferred for the lines are 
not consistent with each other within the formal uncertainties 
(root-mean-square difference of 2.1 kms -1 ; x? c d = 3.5), indi- 
cating that there are small differences. We do not believe that 
this could be due to differences in the behaviour of the two 
Si IV lines, which, as mentioned, arise from the same multi- 
plet. Instead, it might be due to differences in other lines in the 
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Fig. 9. A typical cross-correlation function and the fit applied to the area near the maximum for Si IV 4088.85 A (upper left), 
N II 5676/5679 A blend (upper right), Al III 5696.47 A (lower left) and the interstellar band (lower right). The offset is measured 
in kms . 



two regions; e.g., the 41 16 A region is influenced by the wing 
of the He 14121 Aline. 

Indeed, the velocities inferred from the N III region show 
larger differences with those inferred from the Si IV lines (rms 
of 4.6 kms -1 ; Xrcd = 10)> suggesting that the N III line pro- 
file varies in a somewhat different manner than the Si IV lines. 
This conclusion remains even if one removes the five most dis- 
crepant points, which occur near phase 0.5, when the N III ve- 
locity may be biased towards positive values by the blue-shifted 
absorption feature of HS. 

The radial velocities inferred from the various regions in 
the red spectra are much more discrepant. Again, the velocities 
inferred from regions dominated by lines from the same mul- 
tiplet are more similar to each other (for the two N II regions: 
rms of 9 kms -1 ; xjr cd = 5.4), while the velocities inferred 
from the Al III line are more different (rms of 9.7 kms -1 ; 

XL = ID- 

In summary, we find that for the blue regions, the veloci- 
ties inferred from the different lines track each other closely, 
but not exactly, while for the red region the correlation is much 
less tight. We stress, though, that no systematic differences are 
present. Indeed, below we will find that for all regions con- 
sistent radial-velocity amplitudes are inferred. Also, the differ- 
ences from the best-fit radial-velocity curve will turn out to be 



substantially larger than the differences between the different 
lines, reflecting the fact that, to first order, the lines do behave 
similarly. 

5. The radial-velocity curve of HD 77581 

Ideally, the radial velocities reflect only the motion of the su- 
pergiant companion around the system's center of gravity, and 
are described by a function characterizing a Keplerian orbit 
with known period P, eccentricity e, periastron angle ui, radial- 
velocity amplitude K opt and systemic velocity 7. We fitted the 
radial velocities derived for the different lines, leaving as free 
parameters only K opt and 7Q, since the other parameters are 
const rained very accurately by X-ray timing (Bildsten et al. 
1997| see Tableg). 



In practice, we are confronted with the fact that the indi- 
vidual radial-velocity measurements deviate significantly from 
the best-fit Keplerian orbit. This problem was already recog- 
nized in previous studies of Vela X-l. In Paper II, we found 
that these velocity "excursions" are correlated on a time scale 
of a day and must be non-orbital. This is the reason why we 
chose for the observational strategy to obtain one spectrum per 



Since cross-correlation yields only relative velocities, the fitted 
value of 7 does not reflect the true systemic velocity of the binary. 
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Fig. 11. Comparison of the radial velocities derived from different spectral lines. The difference in radial velocity is shown with 
respect to the radial velocity derived for one of the lines. The label at the top of each figure indicates which lines are involved. 
The left and right figures relate to the blue and red spectral domain, respectively. 




Fig. 10. The radial velocity of the diffuse interstellar band 
(DIB) as a function of orbital phase (filled circles, shifted by 
—75 kms -1 ). For comparison, the radial velocities are shown 
for the N II 5686 A (triangles) and Al III (open circles) re- 
gions; these lines have comparable depth in the red spectra. As 
expected, the velocity of the DIB is consistent with a constant 
value. 



night over many orbits of the system, with the expectation that 
the velocity excursions would average out, thus allowing an ac- 
curate determination of the radial-velocity amplitude. 

we show the radial velocities derived from the 



In Fig. 
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two Si IV and the N III regions separately, as well as the 
radial-velocity curve derived from all these regions combined. 
Overlaid are the best-fit radial-velocity curves; the fitted radial- 
velocity amplitudes are listed in Table || One sees that there 
are large random velocity deviations from the Keplerian orbit, 
with root-mean-square deviations of ~ 7 kms -1 . Superposed 
on these random velocity excursions, there appear to be sys- 
tematic deviations with orbital phase: the velocities are mainly 
above the best-fit curve between phases 0.15 and 0.25 and be- 
tween 0.7 and 0.9, and mainly below the fit between phases 
0.45 and 0.65. We will return to this in Sect. fTi. 



E 




Fig. 12. Radial- velocity curves derived from the two Si IV and 
the N III regions separately, and one derived from the three 
regions combined (labelled "All lines"). The individual radial- 
velocity measurements (filled circles) with their formal error 
bars are shown together with the orbital fit to the curves (drawn 
line; free parameters for the fit are the relative 7-velocity and 
the radial-velocity amplitude K; see Table ||). The dotted line 
reflects the best-fit curve of Fig. 15 
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Fig. 13. Radial-velocity curve derived from the H5 line. Note 
the large deviation towards negative velocities around orbital 
phase 0.5 due to the presence of a photo-ionisation wake. 




Fig. 14. Radial-velocity curves derived from each region sepa- 
rately and from all regions together for the red spectra. 



The radial- velocity curve derived from the HS line (Fig. [l3|) 
shows much stronger systematic deviations with orbital phase, 
in particular an on average large negative velocity excursion 
between phases 0.45 and 0.65 (when the neutron star passes in 



front of the supergiant). As we discussed in Sect. 3.1, the HS 
profile is strongly distorted at these orbital phases, due to an 
additional blue absorption component. This component moves 
towards more negative velocities, also distorting the N III line 
profile. This likely is also the reason for the slight difference 
between the N III and Si IV radial-velocity curves. We did not 



include the region including HS in the final determination of 
the amplitude of the radial- velocity curve. 

In Fig. |l4], we show the radial-velocity curves derived from 
cross-correlation of the N II 5676, 5679 A, the N II 5686 A, 
and the Al III regions, as well as of all regions combined. Fitted 
amplitudes are listed in Table ^. As mentioned, we did not use 
the C III 5710.77 A region because of possible contamination 
with the diffuse interstellar band. 

From the figure, the scatter around the fitted radial-velocity 
curves is much larger than in the blue, with root-mean-square 
deviations of up to 20 krns" 1 (Table |]). On top of this large 
scatter, systematic trends are harder to see, but in general the 
behaviour seems consistent with that in the blue. 

We will return to the effects of the possible systematic ef- 
fects with orbital phase below, after discussing the radial ve- 
locities inferred from ultraviolet spectra, about which there has 
been some controversy, and including these and other velocities 
in a combined velocity curve. 

6. Reanalysis of the IUE spectra 

A large number of ultraviolet spectra of HD 77581 have been 
obtained with the International Ultraviolet Explorer (IUE). In 
Paper II, we used these to determine radial velocities as well. 
We found relatively large uncertainties, of about 2.5 kms -1 , 
but given the large amplitude of the intrinsic variability, they 
nevertheless were a useful supplement to our optical observa- 
tions. 

After our study, a separate analysis of the IUE spectra was 



performed by Stickland et al. (1997). These authors included 



a larger set of spectra and determined velocities by cross- 
correlation with spectra of k Ori rather than with spectra of 
HD 77581 itself. They also used a somewhat different range 
of wavelength regions for calculating the cross-correlation pro- 
files. The radial-velocity amplitude they derived is significantly 
smaller than that found by us; it is also inconsistent with that 
derived from optical spectra in Papers I & II. Clearly, in order 
to have confidence in our results, it is necessary to understand 
the origin of this discrepancy. 

We have contacted Drs. Stickland and Lloyd and they 
kindly provided us with the IUE spectra (reduced using their 
procedures), as well as the procedures and wavelength regions 
used for normalization and cross-correlation. We found that the 
main cause for the discrepancy was that in the Stickland et al. 
cross-correlation procedure for determining the stellar velocity, 
regions around interstellar and wind lines were zeroed in both 
programme-star and comparison-star spectra. Combined with 
slight inadequancies in normalization, this led to the introduc- 
tion of artificial lines at fixed wavelengths in both spectra. In 
consequence, the radial velocities are biased systematically to- 
wards zero and the radial-velocity amplitude is underestimated. 
Dr. Stickland confirmed this conclusion, and reports that after 
correction of the procedure, the radial-velocity amplitude in- 
creases to 22 km s _1 (compared with the 17.8 km s _1 reported 
in their publication), thus solving the discrepancy. 

In order to provide what we hope will be definitive radial 
velocities for the IUE spectra, we decided to re-do the complete 
analysis. For this purpose, we retrieved all 52 SWP spectra of 
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HD 77581 from the archive. These archive spectra have been 
processed using NEWSIPS, the final IUE reduction pipeline. 
We determined radial-velocity differences between all pairs of 
images in the same way as done in Paper II, using the same 
wavelength regions (Table 4 in Paper II), and used these ve- 
locity differences to infer velocities for the individual spec- 
tra. Furthermore, we determined radial velocities relative to the 



Table 3. Radial-velocity amplitudes for the different datasets. 



spectra of k Ori used by Stickland et al. ( |1997[ ): SWP 7745 
(small aperture) and SWP 13776 (large aperture). 

We found that nine spectra of HD 77581 were of insuffi- 
cient quality. For eight of these, the background was abnor- 
mally high (SWP 19041, 19062, 25762, 25763, 25847, 25848, 
25849, 25851), while for the remaining one, the integration was 
very short (1 ks; SWP 25850). We have not used any of these 
further. 

In our analysis, it became clear that many of the cross- 
correlation profi les sho w small excess peaks. These peaks, first 
found by Evans ( 1988[ ) and discussed at some length in Paper II 
(see Fig. 4 in that paper) are thought to result from the presence 
of so-called fixed-pattern noise in the detector. Apparently, the 
new calibration employed in NEWSIPS has not resolved this. 
They occur only if the spectra are taken through the same aper- 
ture. As in Paper II, we corrected for such excess peaks by in- 
cluding an extra Gaussian peak in the fit to the discrete cross- 
correlation profile (if spectra were taken through the same aper- 
ture). 



In Table A. 3 in the Appendix, we list the radial veloc- 
ities derived from the cross-correlations of all spectra with 
each other. We verified that these are consistent with those 
derived from cross-correlation with k Ori. We prefer not to 
use the latter, however, because n Ori has a smaller v sin i 
than HD 77581, and thus is not a good template. We note 



that Stickland et al. (1997) argued that this actually makes it 
a better choice as a template. We do not agree: as one uses a 
template with narrower and narrower lines, approaching delta 
functions, the cross-correlation procedure reduces more and 
more to a weighted sampling of the line profile. Hence, the 
position of the peak of the cross-correlation profile will be- 
come a measurement simply of the position of the deepest 
part of the line profile. In contrast, when one uses a template 
that matches in line width, all information in the line profile 
is used (in a way closely corresponding to % 2 minimalisa- 
tion). We stress, however, that our result does n ot d epend on 
this choice: from the radial velocities in Table A. 3, we find 



K opt = 21.0 ± 1.6 kms- 1 (Table |), while from the veloc- 
ities determined by cross-correlation with the small and large 
aperture spectra of k Ori, we derive K opt = 22.5 ± 2.0 and 
22.9 ± 2.0 kms -1 , respectively. Similarly, our results would 
change very little if we chose the wavelength regions used by 
Stickland et al. (1997 ) instead of the ones from Paper II. 



7. The combined radial-velocity curve and 
systematic effects with orbital phase 

To obtain the best measure of the radial-velocity amplitude, we 
combined the radial-velocity measurements of HD 7758 1 from 
our blue dataset with those derived from the IUE analysis and 
those obtained in Paper II from 13 photographic spectra and 40 



Dataset 


iVd/iVn" 


rms b 


Ampl. c 


K d 


Paper I 


9,11 








21.8 ± 1.2 


Paper II, all data c 


70/ 


JU 


6.5 


9.0 


20.8 ± 1.3 


CCD 




1 1 

1 1 


4.5 


11.3 


21.9 ± 3.1 


Photographic 


13/ 


13 


4.0 


6.8 


21.5 ± 1.9 


Blue, all lines 


77/ 


75 


7.0 


8.2 


22.0 ± 1.0 


Si IV 4089A 






7.1 


8.3 


22.0 ± 1.0 


N III 4097A 






7.4 


8.8 


21.1 ± 1.1 


Si IV4116A 






7.2 


8.5 


22.5 ± 1.0 


Red, all lines 


104/ 


97 


10.1 


15.0 


23.1 ± 1.5 


N II 5676, 5679A 






10.1 


14.8 


23.6 ± 1.5 


N II 5686A 






11.3 


18.2 


19.5 ± 1.9 


Al in 5697A 






11.8 


19.9 


21.2 ± 2.0 


IUE 


43/ 


40 


9.0 


6.4 


21.0 ± 1.6 


CCD, Phot., Blue, IUE 


173/139 


8.4 


8.3 


21.7 ±0.8 



a Total number of spectra and number of nights. 

b Root-mean-square deviation (in kms -1 ) from a Keplerian curve, 

including both measurement uncertainties and intrinsic deviations. 

c Standard deviation (in kms -1 ) of the normally distributed deviation 

amplitudes required to reproduce the rms and used to estimate the 

uncertainty in the radial-velocity amplitude (see text). 

d Inferred radial velocity amplitude (in km s _1 ). 

c Including the velocities derived in Paper II for 26 IUE spectra. These 

are a subset of those analysed in §n. 



CCD echelle spectra. We do not use the velocities derived from 
our red spectra, given their much larger scatter. The resulting 
radial-velocity curve is shown in Fig. O. The best-fit radial- 
velocity amplitude is K op t = 21.7 ± 0.8 kms -1 . In Table pi 
we compare this with the values inferred in earlier sections 
from subsets of the data, as well as with those found previ- 
ously. Clearly, the determinations are all consistent with each 
other. 

7. 1. Search for systematic effects with orbital phase 

Our approach of obtaining radial velocities over an extended 
period of time was motivated by the hope that for a given orbital 
phase the large and highly significant radial-velocity deviations 
from a Keplerian orbit would average out. In other words, that 
just by adding measurements with random deviations, the ac- 
curacy of the derived radial-velocity amplitude, and thus the 
mass of the neutron star, would be improved. If there were to 
be systematic deviations with orbital phase, however, the mea- 
sured mean radial velocity at a given orbital phase would not 
necessarily be equal to the Keplerian radial velocity of the star 
in its orbit. 

We can test the above assumption by verifying whether the 
residuals binned over certain fractions of orbital phase are con- 
sistent with zero. In the lower panel of Fig. O, we show the 
mean of the residuals determined for 9 phase bins, as well as 
the associated errors in the mean, calculated from the standard 
deviation in the phase bins. It is clear that these are not consis- 
tent with zero. The root-mean-square deviation of these phase- 
binned residuals is 2.7 kms" 1 , while only 1.35 kms -1 is ex- 
pected based on the Monte-Carlo simulations. From those sim- 
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Fig. 15. Radial-velocity curve derived from the combined data of this paper and Paper II (top panel). Overdrawn is the Keplerian 
curve that best fits the nightly averages of the data (solid line), as well as the curve expected if the neutron star has a mass of 
1.4 M Q (dotted line; K op t = 17.5 kms" 1 ). The residuals to the best-fit are shown in the middle panel. For clarity, the error bars 
have been omitted. Points taken within one night are connected with lines. In the bottom panel, the residuals averaged in 9 phase 
bins are shown. The horizontal error bars indicate the size of the phase bins, and the vertical ones the error in the mean. The 
dotted line indicates the residuals expected for a 1.4 M Q neutron star. 



ulations, the likelihood that such a high value occurs by chance, 
under the assumption of no orbital-phase related systematic ef- 
fects, is smaller than 0.1%. 

Another piece of evidence that systematic, phase-locked ef- 
fects occur, is that the eccentricity one infers from the optical 
data is 0.18 (and w = 355°). From our simulations, the proba- 
bility to find a value this much larger than the BATSE value is 
less than 1%. 

Understanding the physical origin of these systematic devi- 
ations might allow one to account for them and to further im- 
prove the accuracy of the determination of the radial-velocity 
amplitude. We discuss three possible mechanisms below. 

7.2. Wind effects and evidence for the existence of 
a photoionisation wake 

In the orbital phase range 0.45-0.65, strong velocity excursions 
are observed in the radial-velocity curve obtained for the HS 
line (Sect. In this phase interval, the X-ray source passes 
through the line of sight of the B supergiant. Inspection of the 
HS line shows that at these orbital phases the absorption profile 



includes an additional blue-shifted absorption component, sim- 
ilar to what is observed in time series of spectr al line s formed in 
the supergiant's stellar wind (e.g., Kaper et al. 1994 ). In conse- 
quence, the cross-correlation procedure yields a more negative 
radial velocity for the HS line. 

The origin of this additional blue-shifted absorption com- 
ponent is likely a so-called photo-ionization wake. The X- 
ray photons emitted by the compact X-ray source fully ionize 
the surrounding wind regions, creating an extended Stromgren 
zone in the stellar wind. The presence of such a zone was pre- 
dicted by Hatchett & McCray (1977), and can explain the ob- 
served strong orbital modulation of the ultraviolet resonance 
lines of HD 77581 (e.g., Kaper et al. |1993[). Due to the high 



level of ionization of the plasma contained in the Stromgren 
zone, inside this zone the radiative acceleration of the stel- 
lar wind is quenched, leading to low wind velocities within 
the Stromgren zone. Therefore, at the trailing border of the 
Stromgren zone a strong shock is formed where the relatively 
fast ambient flow meets the slow wind that moved through the 
Stromgren zone. The formation of such a "photo-ionization 
wake" is clearly seen in hydrodynamical simulations (Blondin 
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Fig. 16. Sequence of Ha profiles (left) obtained in February 1996, together with the simultaneously taken red spectra. 



et al. 199C). Additional evidence for the existence of such 



a structure was found in the X-ray light curve of Vela X- 1 
(Feldmeier et al. |1996|). 



On two occassions, in February and May 1996, we obtained 
spectra of strong optical lines formed in the stellar wind to- 
gether with some of the (red) spectra that have been used to 
measure the star's radial velocity. In the period 10-18 February 
1 996 we monitored the Ha line of HD 77581; Fig. |l6| shows the 
dramatic (wind) variability observed in this line. At phase 0.54 
a strong, blue-shifted (—200 to —300 kms -1 ) absorption com- 
ponent appears, which at later phases migrates towards more 
negative velocities. Another blue-shifted absorption compo- 
nent, most probably the remnant of a previous passage of the X- 
ray source, is seen in the first spectra obtained (<f) = 0.09, 0.20, 
0.32). The emission component of the Ha profile is variable 
as well; notice the blue-shifted emission component present at 
phase 0.3-0.4, most likely caused by emission from the dense 
photo-ionization wake when the X-ray source is approaching 
the line of sight. For comparison, Fig. |l^ also displays the red 
spectra obtained within the same nights, where zero velocity 
corresponds to the rest wavelength (in the heliocentric frame) 
of the N II line at 5679.56 A. This line profile clearly demon- 
strates intrinsic variability which cannot be solely due to or- 
bital motion. However, apart from the enhanced blue-shifted 
absorption at (j> = 0.54, there does not seem to be an obvious 
correlation between the variability in the Ha and the N II lines. 
Remarkable is the shallow (incipient emission?) and broad (ad- 



ditional blue-shifted absorption component?) absorption pro- 
file observed at phase 0.2. 

From May 10 to 16, we monitored the He I 5876 A line, 
which is not as strong as the Ha line, but is still sensitive 
to changes in the base of the stellar wind. The systematic 
changes (Fig. |l7|) strongly resemble those observed in the H/3 
line (Kaper et al. 1994) and are similar to the variations found 



in Ha. Comparison with the red N II spectra shows that at 
= 0.34 both the He I and the N II (and Al III) lines are shal- 
lower, while at cj) = 0.57 all lines are deeper. The blue-shifted 
absorption is enhanced at phases 0.46 and 0.56. 

We conclude that the systematic deviation in radial velocity 
of the HS line in the phase interval 0.45 — 0.65 is caused by 
the presence of a photo-ionization wake. Weaker spectral lines 
are much less strongly, though still measurably affected in this 
phase interval; stronger lines like Ha show large deviations in 
line shape over a much wider range of the orbit. We decided to 
leave out the radial-velocity measurements based on the HS line 



in our further analysis. In Sect. 7.5, we investigate the impact 



of the points between phase 0.45 and 0.65 on the determination 
of the radial-velocity amplitude. It turns out to be small. 

7.3. Tidal distortion 

Tidal deformation could also cause systematic deviations of the 
radial-velocity curve, since for a deformed star 'center-of-light' 
radial-velocity measurements do not necessarily reflect the ac- 
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Fig. 17. Sequence of He I profiles (left) obtained in May 1996 compared to the simultaneously taken red spectra. 



tual center-of-mass velocity. To investigate the effect this might 7. 4. Non-radial pulsations 



have on inferred orbital parameters, Van Paradijs et al. (1977b) 



made model calculations for a star with properties appropri- 
ate for HD 77581. For simplicity, they assumed a circular orbit 
and synchronous rotation of the primary, so that the shape of 
the star was given by a single Roche equipotential surface at 
all phases. Using theoretically predicted line profiles, includ- 
ing the Si IV 4089 line, they find that appreciable systematic 
effects can occur, differing from ion to ion, with the main ex- 
pected deviations being a positive one around orbital phase 0. 1 
and a negative one around phase 0.9. As a result, an apparent 
eccentricity may arise and the radial-velocity amplitude may 
be overestimated, in unfavourable cases by up to 30%. 



The systematic deviations we encounter in the radial- 
velocity curve behave differently. The larger deviations occur 
around orbital phases 0.2 and 0.8, and at both phases the de- 
viations are directed towards positive velocities. Thus, it seems 
the effects of the tidal deformation are not as large as predicted. 
One should keep in mi nd, how ever that the assumptions made 
by Van Paradijs et al. ( 1977b ) do not correspond to the actual 
situation in Vela X-l: the orbit is eccentric and HD 77581 ro- 
tates sub-synchronously (by a factor of about 2/3, Zuiderwijk 



1995). Given the possibly large effects on the inferred parame- 



ters, it would be worthwile to reconsider these issues. 



Apart from distorting the shape of the star, tidal forces can 
excite (non-radial) pulsations. With a neutron star in a close 
eccentric orbit (the distance between the neutron star and the 
non-synchronously rotating surface of HD 77581 is about half 
a stellar radius) one can expect that tidal waves ar e exci ted at 
the surface of HD 77581 (e.g., Witte & Savonije |1999| ). The 
neutron star's orbital frequency may resonate with certain non- 
radial modes, perhaps phase-locked with the orbit. The proper- 
ties of these modes depend on the details of the internal struc- 
ture of the supergiant, which has undergone a complex history 
of binary evolution, including a phase of mass transfer and a 
nearby supernova explosion forming the neutron-star compan- 
ion. 

Radial and low-degree non-radial pulsations can change 
the shape of photospheric absorption line profiles significantly 
(e.g., Vogt & Penrod 1983). This is due to local temperature 



variations induced by the pulsations and/or local Doppler shifts 
related to the 3-dimensional pulsation velocity field. 

The question is whether such tidally induced (or self- 
excited) pulsations, when present in HD 77581, are detectable 
through inspection of the detailed shape of the line profiles 
and, if so, could result in systematic deviations injhe radial- 
velocity curve. According to Telting & Schrijvers (1997), low- 
order pulsation modes (with degree £ < 3) can lead to de- 
tectable radial-velocity variations and changes in the shape of 
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the line-profiles. For instance, moment analyses of the spec- 
tral line profiles (Balona 1986, Aerts et al. 1992| ) have been 
successfully used to detect low-order modes in single pulsators 
such as the j3 Cephei stars. 

If the presence of non-radial pulsations is the physical cause 
for the observed deviations in radial velocity, knowledge of the 
pulsation mode(s) might provide a tool to "correct" the mea- 
sured radial velocities and to improve the accuracy of the de- 
termination of the Keplerian orbit of HD 77581. However, a 
mode can only be identified if the time sampling of the data is 
sufficient to resolve its periodicity (in practice several spectra 
per hour and continuous - 24 hours per day - coverage). The 
time sampling of the spectra used in this paper (one spectrum 
per day) is poor; it is sufficient to cover the change in radial 
velocity due to orbital motion, but likely insufficient to identify 
pulsation modes. Another disadvantage of the current dataset 
is the limited spectral coverage; this excludes the possibility 
to combine many different spectral lines in order to increase 
the signal-to-noise ratio and to separate line changes induced 
by the temperature- and velocity distribution associated with a 
given pulsation. 

A preliminary moment analysis performed on the line pro- 
files comprised in our spectra did not yield a useful result. This 
is an indication that any putatuve non-radial pulsation modes 
are not phase-locked to the orbital motion and/or do not persist 
for a long period of time. A dataset obtained with better time 
sampling and spectral coverage is required to further investi- 
gate the occurrence of non-radial pulsations in the atmosphere 
of HD 77581. 



7.5. Orbital-phase intervals excluded 

In the H<5 profile, we found clear evidence for an ionisation 
wake in the phase interval between 0.45 and 0.65. It may be 
that the photoionisation wake has a (less obvious) effect on 
the other spectral lines as well. We can verify to what extent 
the exclusion of this phase interval would affect the amplitude 
of the radial-velocity curve (and thereby the measurement of 
the neutron-star mass). Similarly, in previous studies it was 
found that the largest deviations occurred near velocity mini- 
mum (Paper II and references therein). From Fig. |l5[ one sees 
that the phase bin at velocity minimum again has one of the 
largest and most significant deviations from zero. 

In order to see what is the effect of these phase intervals on 
the solution, we decided to exclude, like in Paper II, all possi- 
ble one ninth and two ninth cycle wide orbital-phase intervals 
(i.e., approximately one and two-day wide) from the fit to the 
radial-velocity curves, starting from orbital phase zero until one 
with a step of one ninth. Subsequently, we fitted the remaining 
points and derived the amplitude of the radial-velocity curve, as 
well as the uncertainty estimated using our Monte-Carlo sim- 
ulations. The resulting amplitude is shown as a function of the 
excluded phase intervals in Fig. |l8[ Also shown is the rms of 
the phase-binned residuals. If a specific phase interval were re- 
sponsible for virtually all of the phase-locked deviation, one 
would expect that the binned residuals would become consis- 
tent with zero when that phase bin is excluded. From the fig- 
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Fig. 18. Results for fits to the radial-velocity curve with phase 
intervals of width 1/9 (left) and 2/9 (right) excluded. The points 
in the top panel show the best-fit radial-velocity amplitudes 
found when the phase range indicated by the horizontal error 
bars is excluded. The short-dashed contours indicate the 95% 
and 99% confidence range of the variation expected from our 
simulations (for which it is assumed that the variability of the 
star is independent of orbital phase.) The long-dashed lines in- 
dicate the best-fit amplitude to all data and the associated 95% 
confidence uncertainties. The panels in the second row show 
the root-mean-square deviation of the phase-binned residuals. 
The short-dashed lines indicate the levels below which 50%, 
95%, and 99% of the simulated data are contained. In the bot- 
tom two rows of panels, the best-fit values for the eccentricity 
and periastron angle are shown. The long-dashed lines indicate 
the accurate values inferred from the BATSE measurements 
(Table Q), the short-dashed contours the 95% and 99% ranges 
expected from the simulations. 
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ure, one sees that excluding the bins near inferior conjunction 
of the neutron star, when the ionisation wake is most important, 
and near velocity minimum have the largest effect on the rms. 
Excluding inferior conjunction, the inferred radial-velocity am- 
plitude is not affected much, since the velocities are at a zero 
crossing, while excluding velocity minimum, a much larger ve- 
locity amplitude is inferred; indeed, excluding this bin leads to 
the only change in radial-velocity amplitude whose probability 
of occurring by chance, as inferred from the simulations, is less 
than 1%. 

7.6. Fitting eccentricity and periastron angle 

In another attempt to determine whether any specific bin is re- 
sponsible for the phase-locked systematic effects, we repeated 
the above simulations but with this time the eccentricity and 
the periastron angle as free parameters in the fitting procedure. 
In principle, of course, one should reproduce the values found 
from the X-ray analysis. The bottom panels of Fig. [l8| show the 
resulting eccentricity and periastron angle. We find that when 
early phase intervals are excluded, the inferred eccentricity re- 
mai ns far larger than the true value, measured by Bildsten et 
al. ( 1997 ; Table Q). When phase intervals around velocity min- 
imum are excluded, however, the inferred eccentricity becomes 
much closer to the true value. 



7. 7. The uncertainty in the radial-velocity amplitude 

The above results suggest that the main systematic effects on 
the radial-velocity curve occur near inferior conjunction and 
at velocity minimum. If so, the radial-velocity amplitude in- 
ferred from all the data likely is an underestimate of the true 
value. However, while we believe the velocity deviation at in- 
ferior conjunction can be understood - in terms of a photo- 
ionisation wake - we do not have a good idea for the cause of 
the other systematic deviations. Lacking this, it is difficult to 
estimate the uncertainty that is introduced. Probably best is to 
assume there is no bias to high or low radial-velocity ampli- 
tude, and base the estimate on the increase in the uncertainties 
required to match the observed root-mean-square deviations of 
the phase-binned residuals. This is about a factor two, and cor- 
responds to a similar increase in the uncertainty on the final 
radial-velocity amplitude. In consequence, our final estimate is 
not much more precise than values given previously, in which 
the influence of systematic effects with orbital phase was not 
taken into account. We conclude that our best present estimate 



Combined with the X-ray mass function 
Ml sin 3 i An 2 



f: 



(axsiniy = 19.74 ± 0.07 M 



is K, 



opt 



21.7 ±1.6 km s _1 . 



8. The mass of Vela X-1 

From the ratio of the radial-velocity amplitudes of the optical 
and X-ray components, we infer a mass ratio 

q = M x /M opt = Kopt/Kx = 0.078 ± 0.006, 

after inserting the value of the radial-velocity amplitude of the 
X-ray source (derived from X-ray pulse-timing analysis) 

K x = (27r/P)a x sinz(l - e 2 )~ 1/2 = 278.1 ± 0.3 km s" 1 . 



{1 + l/q) 2 GP 2 
we infer 

Mx sin 3 i = 1.78 ± 0.15 M Q , and 
M opt sin 3 i = 22.9 ± 0.3 M Q . 

The above numbers are lower limits to the actual masses, 
since the inclination can only be smaller than 90°. As in 
Paper II, we use the method introduced by Joss & Rappaport 
(1983) to estimate the inclination i, the semi-major axis a, 
the radius i? pt and mass M opt of the companion, the co- 
rotation factor f co , and the neutron-star mass Mx, from the 
observed values of the orbital period P, eccentricity e, peri- 
astron angle u, projected semi-major axis of the neutron star 
ax sin i, eclipse duration 9 CC \, radial-velocity amplitude K opt , 
and projected rotational velocity v mt sini (see Table The 
corresponding uncertainties are found by means of a Monte- 
Carlo error propagation technique, in which a large number of 
trial evaluations are made for parameters drawn from the ob- 
served uncertainty distribution (which are assumed to be nor- 
mal, except for 9 CC \, for which a uni form distribution is used). 
Following Joss & Rappaport ( 1983 ), we assume that the star 
is close to filling its Roche lobe, with a volume filling factor 
between 90 and 100% at periastron. 

The results are listed in Table for the mass of the neutron 
star, we find 



M x = 1.86 ± 0.16 M Q . 

Here we should add a note of caution: as in Paper II, for a 
substantial fraction (37%) of the trials, the randomly-chosen 
parameters are inconsistent with each other in the sense that 
the eclipse duration cannot be reproduced for any inclination. 
We discard these trials, but note that as a result, the inferred 
parameters are based on effective distributions of the input pa- 
rameters which no longer reflect the actual ones. Most likely, 
this is because the assumption that the star is in instantaneous 
hydrostatic equilibrium is incorrect (see Paper II for a more ex- 
tensive discussion). 

Given the above uncertainty, and given that for constraints 
on the neutron-star equation of state, the most interesting result 
is a firm lower limit to the mass of the neutron star, we prefer 
to give lower limits in terms of Mx sin 3 i. For these, we use the 
one-sided 95% and 99% confidence lower limits to the radial- 
velocity amplitude, of 19.1 and 17.9 kms , respectively, for 
which we infer M x sin 3 i > 1.54 and 1.43 M Q , respectively. 

9. Ramifications 

The best estimate of the mass of Vela X-l is 1.86 M Q . 
Unfortunately, no firm constraints on the equation of state 
are possible, since systematic deviations in the radial-velocity 
curve do not allow us to exclude a mass around 1.4 Mq as 
found for other neutron stars. Therefore, we will not discuss 
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Table 4. Parameters of the Vela X-l/HD 77581 system. 



Parameter 



Value a 



Unit 



Ref. 



Observed (la uncertainties) 



ir/2 
-forb 

Ox sin i 
e 

Vrot sin i 



K, 



opt 



2448895.7186(12) 

8.964368(40) 

113.89(13) 

0.0898(12) 

152.59(92) 

30-36 

116(6) 

21.7(16) 



JD 

day 
lt-s 



kms 1 
kms -1 



Bildsten et al. 



See Paper I 



1997 



Zuiderwijk 
This paper 



1995 



Derived (la 


uncertainties) 




M x sin 3 i 


1.78 ±0.15 


M 


M opt sin 3 i 


22.9 ±0.3 


M 


Inferred (95% confidence ranges' 1 ) 




fco 


69 +0 ' 08 




i 


>73 


o 


a 


53.6t£;J 


R Q 


-Ropt 


30.41J;? 


Rq 


M opt 


23.8+.^ 


M e 


M x 


1.86l° ;g 





a Time of mean longitude 90° . 
b 99% confidence range. 

c For the inferred quanties, 95% confidence ranges are listed, since 
the uncertainties on many parameters have highly non-normal dis- 
tributions. The one exception is the mass of the neutron star, whose 
uncertainty distribution is close to normal. Hence, the associated la 
error is about 0.16 Mq. 



the equation of state here, but only address evolutionary ram- 
ifications of what still most likely is a massive neutron star, 
assuming that the equation of state is sufficiently stiff for such 
a massive neutron star to exist. 

One reason that the high mass is peculiar (if true), is that the 
most accurately determined neutron-star masses fall in a sur- 
prisingly small range: for the binary radio pulsars with neutron- 
star companions, the masses are in the range 1.338-1.442 Mq. 
Indeed, Thorsett & Chakrabarty (1999) found that all binary 
radio pulsars have masses which are consistent with a nor- 
mal distribution around 1.35 Mq, with the remarkably small 
spread of 0.04 Mq. Furthermore, for all X-ray binaries other 
than Vela X-l for which dynamical masses are available, the 
neutron-star masses are consistent with the same value (Van 



Kerkwijk et al. 1995c ) 



Obviously, since for no equation of state masses substan- 
tially lower than 1.4 Mq are excluded, at least the lower limit 
to the above narrow range must be set by the way in which the 
neutron star is formed, i.e., by the physics of supernova explo- 
sions and the evolution of stars massive enough to reach core 
collapse. This leads one to wonder whether it could not be the 
formation mechanism on its own that leads to the narrow range 
in mass, just as for white dwarfs the formation leads to masses 
mostly within a very narrow range around 0.6 Mq, well below 
the Chandrasekhar mass. 



If the formation indeed leads to a narrow range of masses, 
high masses could arise only if substantial amounts of mat- 
ter are accreted. This is not expected to have occurred for 
any of the radio or X-ray pulsars with accurate mass esti- 
mates. For low-mass X-ray binaries, in which substantial ac- 
cretion is expected to have occurred, masses around 2 Mq 
have been suggested (see, e.g., Orosz & Kuulkers [1999] for 



an analysis of Cyg X-2, and Zhang et al. [ 1997] for inferences 
based on quasi-periodic oscillations). All these estimates, how- 
ever, rely to greater or lesser extent on unproven assumptions. 
Furthermore, for the putative descendants, radio pulsars with 
white dwarf companions, there is no evidence for such high 
masses (Thorsett & Chakrabarty 1999 and references therein). 
From recent evolutionary calculations, Timmes et al. 



( 199C ) indeed expect that massive stars in close binaries, which 
explode as Type lb supernovae, give rise to initial neutron star 
masses in a narrow range around 1.3 Mq. This values does 
not include subsequent mass accretion from a reverse shock 
or from a massive component in a binary system, and Timmes 
et al. expect that the final masses could be somewhat higher. 
Interestingly, for single stars, which explode as Type II super- 
novae, they found a bimodal distribution of initial neutron-star 
masses, with narrow peaks at 1.27 and 1.76 Mq. As mentioned, 
they did not find a bimodal distribution for stars in close bina- 
ries, but at present it is not clear whether this result will hold 
(Woosley 2000, private communication). If stars in close bi- 
naries turn out to be more similar to single stars after all, one 
could assign most neutron stars to the first peak, and Vela X-l 
to the second. 

We end by noting that even if neutron stars are formed with 
masses in two peaks, it might be rather unlikely to find a mas- 
sive neutron star in a close double neutron star binary. This 
is because such systems require a common-envelope stage, in 
which a merger can only be avoided if the initial orbit was very 
wide. Stars massive enough to form a massive neutron star, 
however, likely do not evolve through a red-giant phase, and 
a common-envelope phase would occur only for rather close 
orbits, for which the binary would merge. 
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Appendix A: Log of observations 

The tables in this Appendix listing the log of observations 
of the blue CES spectra (Tab. A.l), the red CES spectra 
(Tab. |A.2| ), and the IUE spectra (Tab. A. 3) are only available 
in electronic form at the CDS via anonymous ftp to cdsarc.u- 
strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi- 
bin/qcat?J/A+A/. To provide an impression of the information 
contained in these tables, only the first five lines of each table 
are displayed here. 
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Table A.l. Observing log and radial velocities for the blue spectra (central wavelength 4102 A). The second digit of the identifica- 
tion (first column) specifies the CCD used (b4 for CCD#34 and b8 for CCD#38). In the second column the barycentric-corrected 
Julian Date of each observation is listed, and in the third the orbital phase calculated using the ephemeris of Bildsten et al. 



(1997; see Table |4|). The following four columns list radial velocity and associated uncertainty (in kms -1 ) for each spectral 
region separately, as derived by means of cross-correlation. The last columns lists radial velocity and uncertainty derived from 
the cross-correlation of the ensemble of the two Si IV and the N III regions. 



Ident 


BJD m id 


4>orb 


Si IV 4089 


N in 4097 


H5 4102 


Si iv 41 16 


Si/N ensemble 


b8001 


2450006.8514 


0.9502 


-5.99±0.82 


-9.60±1.47 


-5.71±0.77 


-3.87±1.05 


-6.27±0.51 


b8002 


2450006.8660 


0.9518 


-5.31±0.83 


-8.05±1.43 


-6.52±1.06 


-3.89±1.04 


-5.63±0.51 


b8004 


2450008.8590 


0.1741 


23.11±0.97 


21.90±1.75 


29.18±0.81 


22.66±1.26 


22.69±0.58 


b8005 


2450009.8532 


0.2850 


18.14±0.71 


14.68±1.36 


24.07±1.15 


15.72±1.04 


16.79±0.45 


b8006 


2450010.8543 


0.3967 


10.76±0.43 


10.87±0.81 


11.75±0.94 


9.90±0.62 


10.23±0.31 



Table A.2. Observing log and radial velocities of the red spectra (central wavelength 5965 A). Columns as in Table 4.1 



Ident 


BJD m id 


4>orb 


N II 5676/9 


N II 5697 


Al m 5686 


C in 5705 


N/Al ensemble 


r8002r 


2450006.8545 


0.9505 


-10.55±1.83 


-7.98±3.97 


-11.61±2.62 


-2.47±4.25 


-10.48±1.38 


r8001r 


2450012.8401 


0.6182 


-8.96±1.82 


-1.77±3.46 


-4.75±2.45 


-1.56±3.78 


-7.65±1.29 


r8003r 


2450013.8483 


0.7307 


-19.48±1.88 


-20.12±3.48 


-18.55±2.93 


-11.60±4.74 


-19.25±1.38 


r8004r 


2450013.8629 


0.7323 


-17.66±1.72 


-17.13±2.95 


-21.05±2.39 


-9.17±3.32 


-17.92±1.22 


r8005r 


2450014.8585 


0.8433 


-7.36±1.59 


-10.61±2.69 


-7.93±1.79 


-1.39±3.77 


-7.43±1.16 



Table A.3. IUE observations of HD 77581. 



SWP 


JDmid. cxp. 


Exp. time 


Orbital 


Velocity b - c 


number 


-2440000 


(min.) 


phase a 


(kms -1 ) 


1442 


3628.708 


180 


0.451 


0.0 ± 1.5 


1488 


3634.118 


150 


0.055 


5.1 ± 1.9 


2087 


3712.994 


125 


0.854 


-28.6 ± 2.0 


2390 


3745.203 


108 


0.447 


3.9 ±2.8 


3499 


3843.592 


120 


0.422 


7.9 ±2.7 



a Using the ephemeris of Bildsten et al. 1997 (see Table ^). 
b Relative to SWP 1442; corrected for shift of interstellar lines. 
c Quoted are la errors. 
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